home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_200 / 247_02 / pollard.c < prev    next >
Text File  |  1989-04-17  |  4KB  |  138 lines
  1. /*
  2.  *   Program to factor big numbers using Pollards (p-1) method.
  3.  *   Works when for some prime divisor p of n, p-1 has only
  4.  *   small factors.
  5.  *   See "Speeding the Pollard and Elliptic Curve Methods"
  6.  *   by Peter Montgomery, Math. Comp. Vol. 48 Jan. 1987 pp243-264
  7.  */
  8.  
  9. #include <stdio.h>
  10. #include "miracl.h"
  11.  
  12. #define LIMIT1 10000        /* must be int, and > MULT/2 */
  13. #define LIMIT2 300000L      /* may be long */
  14. #define MULT   2310         /* must be int, product of small primes 2.3.. */
  15.  
  16. main()
  17. {  /*  factoring program using Pollards (p-1) method */
  18.     int phase,m,iv,pos,btch;
  19.     int *primes;
  20.     long i,p,pa;
  21.     big b,q,n,t,bw,bvw,bd,bp;
  22.     static big bu[1+MULT/2];
  23.     static bool cp[1+MULT/2];
  24.     mirsys(50,MAXBASE);
  25.     b=mirvar(2);
  26.     q=mirvar(0);
  27.     n=mirvar(0);
  28.     t=mirvar(0);
  29.     bw=mirvar(0);
  30.     bvw=mirvar(0);
  31.     bd=mirvar(0);
  32.     bp=mirvar(0);
  33.     primes=gprime(LIMIT1);
  34.     for (m=1;m<=MULT/2;m+=2)
  35.         if (igcd(MULT,m)==1)
  36.         {
  37.             bu[m]=mirvar(0);
  38.             cp[m]=TRUE;
  39.         }
  40.         else cp[m]=FALSE;
  41.     printf("input number to be factored\n");
  42.     cinnum(n,stdin);
  43.     if (prime(n))
  44.     {
  45.         printf("this number is prime!\n");
  46.         exit(0);
  47.     }
  48.     while (gcd(n,b,t)!=1) incr(b,1,b);
  49.     phase=1;
  50.     p=0;
  51.     btch=50;
  52.     i=0;
  53.     printf("phase 1 - trying all primes less than %d\n",LIMIT1);
  54.     printf("prime= %8ld",p);
  55.     forever
  56.     { /* main loop */
  57.         if (phase==1)
  58.         { /* looking for all factors of (p-1) < LIMIT1 */
  59.             p=primes[i];
  60.             if (primes[i+1]==0)
  61.             { /* now change gear */
  62.                 phase=2;
  63.                 printf("\nphase 2 - trying last prime less than %ld\n",LIMIT2);
  64.                 printf("prime= %8ld",p);
  65.                 power(b,8,n,bw);
  66.                 convert(1,t);
  67.                 copy(b,bp);
  68.                 copy(b,bu[1]);
  69.                 for (m=3;m<=MULT/2;m+=2)
  70.                 { /* store bu[m] = b^(m*m) */
  71.                     mad(bw,t,t,n,n,t);
  72.                     mad(bp,t,t,n,n,bp);
  73.                     if (!cp[m]) continue;
  74.                     copy(bp,bu[m]);
  75.                 }
  76.                 power(b,MULT,n,t);
  77.                 power(t,MULT,n,t);
  78.                 mad(t,t,t,n,n,bd);     /* bd = b^(2*MULT*MULT) */
  79.                 iv=p/MULT;
  80.                 if (p%MULT>MULT/2) iv++,p=2*(long)iv*MULT-p;
  81.                 power(t,(2*iv-1),n,bw);
  82.                 power(t,iv,n,bvw);
  83.                 power(bvw,iv,n,bvw);   /* bvw = b^(MULT*MULT*iv*iv) */
  84.                 subtract(bvw,bu[p%MULT],q);
  85.                 btch*=10;
  86.                 i++;
  87.                 continue;
  88.             }
  89.             pa=p;
  90.             while ((LIMIT1/p) > pa) pa*=p;
  91.             power(b,(int)pa,n,b);      /* b = b^pa mod n  */
  92.             decr(b,1,q);
  93.         }
  94.         else
  95.         { /* phase 2 - looking for one last prime factor of (p-1) < LIMIT2 */
  96.             p+=2;
  97.             pos=p%MULT;
  98.             if (pos>MULT/2)
  99.             { /* increment giant step */
  100.                 iv++;
  101.                 p=(long)iv*MULT+1;
  102.                 pos=1;
  103.                 mad(bw,bd,bd,n,n,bw);
  104.                 mad(bvw,bw,bw,n,n,bvw);
  105.             }
  106.             if (!cp[pos]) continue;
  107.             subtract(bvw,bu[pos],t);
  108.             mad(q,t,t,n,n,q);  /* batching gcds q=q.t mod n */
  109.         }
  110.         if (i++%btch==0)
  111.         { /* try for a solution */
  112.             printf("\b\b\b\b\b\b\b\b%8ld",p);
  113.             gcd(q,n,t);
  114.             if (size(t)==1)
  115.             {
  116.                 if (p>LIMIT2) break;
  117.                 else continue;
  118.             }
  119.             if (compare(t,n)==0)
  120.             {
  121.                 printf("\ndegenerate case");
  122.                 break;
  123.             }
  124.             printf("\nfactors are\n");
  125.             if (prime(t)) printf("prime factor     ");
  126.             else          printf("composite factor ");
  127.             cotnum(t,stdout);
  128.             divide(n,t,n);
  129.             if (prime(n)) printf("prime factor     ");
  130.             else          printf("composite factor ");
  131.             cotnum(n,stdout);
  132.             exit(0);
  133.         }
  134.     } 
  135.     printf("\nfailed to factor\n");
  136. }
  137.  
  138.